Probing QCD Phase Structure by Baryon Multiplicity Distribution 
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The canonical partition functions Z n and the number distributions P n which are obervable in 
experiments, are related by a single parameter, the fugacities £ = exp(/i/T). With the charge 
parity invariance, Z n and £ can be determined. Thermodynamic quantities such as the number 
density susceptibility and the kurtosis are then calculated from the grand canonical partition function 
T) = J2 Zn{T)£ n , (n = — iV max , • • • , iV max ), for any chemical potential fi, although the region 
over which the results are reliable for these quantities is constrained by iVmax- We then calculate the 
Lee- Yang zeros, which are the zeros of Z{£) in the complex fugacity plane, as poles of dlog Z(£)/d£ 
by using the Cauchy integral theorem. With the help of a multiple precision library, this method 
provides any precision required without misidentification of the zeros. 

We analyse Z n from the net-proton number distributions recently measured at the Relativistic 
Heavy Ion Collider (RHIC) by assuming the net-proton number is approximately propotional to that 
of the baryon after the freeze-out, and calculate the moments. We also evaluate the Lee- Yang zero 
structures obtained from RHIC data and compare them with those obtained from lattice quantum 
chromodynamics (QCD) calculations. 

Possible regions of QCD phase transition lines are estimated from the thermodynamics quantities 
and the Lee- Yang zeros. We discuss how the limited JV max in both experimental and numerical 
studies affects the reliability of the thermodynamic results and Lee- Yang zeros. 
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When temperature and density are varied, QCD is ex- 
pected to have a rich phase structure [l| ■ One of the most 
important challenges in particle and nuclear physics is to 
experimentaly discover these phases that only appear un- 
der extreme conditions and to theoretically understand 
the nature of these phases. Such achievements would not 
only deepen our understanding of QCD but also extends 
our knowledge of the early universe and compact stars. 

The Relativistic Heavy Ion Collider (RHIC) was built 
to explore properties of QCD matter 0, and has pro- 
vided much valuable information 3, 4]. In addition, many 
phenomenological studies have also focused on the phases 
of QCD matter [5|. Although lattice QCD provids very 
precise information about QCD at finite temperatures [6|, 
thirty years of numerical studies of lattice QCD at finite 
baryon number density are still far from yielding con- 
clusie results because of a notorious sign problem 0, H| . 

Recently net-proton multiplicity measurements at 
RHIC are gaining attention ({J [l(| because they enable 
the QCD phase diagram to be explored [Hj |. In these 
measurements, the colliding energy is varied, and the 
trajectory (T,hb) in the QCD phase diagram may pass 
near the critical region. Even-by-event fluctuations are 
expected to encompass the critical point, where the cor- 
relation length changes. In particlular, conserved quan- 
tities such as charge or baryon number may reveal pos- 
sible correlations that existed inside the system before 
hadronization. (see Ref.[12| and references therein). 

In this paper, we propose a method by which, even 
without experimental data near the phase transition, the 
observed number distribution data can be used to obtain 
more information about the QCD phase transition, and 
even without direct lattice QCD calculations at physical 
chemical potential regions, we can use the data to un- 
derstand experimental multiplicity through the canonical 



partition functions. 

The grand partition function Z and the canonical par- 
tition function Z n are related according to 



Z(Z,T)=Tre-K H -> t *) 



+N ma . 



n=-N„ 



Zn{T)C 



(1) 



where £ = exp(/i/T). Here we assume that the number 
operator N commutes with H, that is, N is a conserved 
quantity. Because of the charge-parity symmetry, Z n sat- 
isfies 



(2) 



The multiplicity distributions P n observed in experi- 
ments are related to Z n as 



Pn(0 = Z n C- 



(3) 



Using Eqs.© and ©, we can determine £ and Z n . 
FigJT] shows the obtained £ togather with that by a 
freeze-out analysis in Ref.|13j. From the grand parti- 
tion functions Eq.([Ij). the moments are evaluated using 

A fc (0^(^) fe iog^(£). 

We assume that the net-proton multiplicity data are 
approximately propotional to those of the baryon. The 
proportionality factor has no effect on any of the results 
reported in this paper. This approximation is justified if 
(i) after the freeze-out, the net-proton number is effec- 
tively constant, or (ii) a created fireball is approximately 
iso-neutral. 

The quantity -/V max in Eq.([T]) is finite because of mea- 
surement statistics (simulation statistics and finite vol- 
ume) in the experiments (lattice QCD). The finite nature 
of JVmax puts an upper bound on the chemical potential 
for which the calculation is valid. To estimate the effect 
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FIG. 1: Fugacity £ = exp(/i/T) as a function of the 
colliding energy y/s. Plotted arc those obtained from 
RHIC experiments (blue crosses) and the freeze-out 
results reported in Ref . fl~3j . 
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FIG. 2: The number susceptibility Aa/iV max , at 
y^s = 27 as a function of fi/T. The upper curve is 
obtained by raising the final three 
Z n , (n = V max - 2, Y max - 1, A max ) by 15%. The lower 
curve results from dropping the final two Z n terms. 



of finite N max , we test two cases: (1) the values of the 
final three Z n (n = N mecK -2,N max -l,N max ) are increased 
by 15% and (2) the final two Z n (n = N max -l,N max ) 
are set to zero. As an example, we plot the number 
susceptibility A 2 /A max in FigJ2] as a function of fi/T at 
y/s = 27 GeV for these two cases together with the one 
constructed from the full Z n . 

Let us suppose that, as /i is increased, we encounter 
a phase transition or a cross-over transition, and cross 
it. In this case, we would expect the structure of the 
moments to be rough in this area. At the peak position 
of the lower curve indicated by an arrow in Figj2] the 
center line continues to increase. We write this value 
as p(T), and any transition may occur for (J, > p. In 
other words, this position p is a candidate for the lower 
bound of the real susceptibility peak. We then check the 
behavior of A3/A2 and A4/A2; Although higher moments 
have large effects because of the finite Nm&x, they may 
work as a good tool to detect the transition region of 




FIG. 3: Ration of the moments A4/A2 as a function of 
(x/T at ^/s =19.6GeV. The upper and lower curves are 
constructed by the same procedure as in Fig|5] 
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FIG. 4: Potential regions of hidden phase transition. 
The regions are estimated as complement of the areas 
where no transition is evident. 



the QCD phase diagram 14} . The former ratio indicates 
no significant structure for fi < p. Around the freeze- 
out points, A4/A2 ~ 1 (Poisson) and it decreases to the 
negative values as /i increases. SeeFigH InRefs.[l!-tTz|, 
it is argued that the negative value of X4/X2 is a signal of 
reaching the phase transition. We take the points where 
the lower curve of A4/A2 becomes negative (indicated by 
an arrow in Figj3]) as an another candidate of the lower 
bound p(T) and show them in FigUJ 

Next we extend the fugacity £ to complex values. Lee- 
Yang zeros (LYZs) ai are zeros of the grand partition 
function Z(£) in the complex fugacity plane; that is, 
Z(£) — rj(£ — on) . The distribution of these zeros re- 
flects the phase structure of the corresponding statistical 
system. Lee and Yang argued that for a finite system, 
no zeros appear on the real axis. In the thermodynamic 
limit, the number of zeros becomes infinite and the zeros 
coalesce onto one-dimensional curves [T^- For the first- 
order (second-order) phase transition, the coalescing ze- 
ros cross (pinch) the real positive £ axis[l9|. whereas for 
the cross-over, they do not reach to the real axis. 

Therefore it is important to reliably determine the 
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FIG. 5: Schematic of cBK contours in the 
divide-and-conquer search for residues. 



LYZs; that is, all zeros must be found without ambiguity, 
and their positions in the complex fugacity plane must 
be determined with high accuracy. We obtain the LYZs 
as follows: We first map the problem to a calculation of 
the residue of Z 1 /Z: 

Z'/Z = (4) 

In Eq.Q, the left hand side is integrated along a contour, 
and the residues inside the contour are summed according 
to Cauchy's theorem. Because of the symmetry o 
— /z, if ai is a solution, i.e., Z{a{) — 0, then 1/ai is also 
a solution. Therefore we only need to search for residues 
inside the unit circle. 

Fig[S] shows cut Baum-Kuchen (cBK) shape contours 
used in the study. Starting from < 8 < 2ir, and 
< r < 1 in polar coordinates, the region is divided into 
four pieces, and the Cauchy integral is evaluated over 
each section. This divide-and-conquer process is iterated 
up to as many times as needed (here 20 times). When no 
residue is found in a section, no further divisions are ap- 
plied to that section. At each divide-and-conquer level, 
we check the conservation of the residue sum. The tech- 
nical details of this approach will be presented elsewhere. 

All calculations were performed using the multiple- 
precision package, FMlib[20|], and the number of signif- 
icant digits was 50 - 100. We checked that residues in- 
volved in the analyses presented herein are 1 + e with 
|e| < 10- 3 . 

With this algorithm, we can construct LYZ diagrams 
from Z n . First, we study the LYZs of the lattice QCD 
simulation. Here we do not distinguish u and d quarks. 
Details of calculating Z n with a lattice QCD simulation 
are given in [2lj |. where the fugacity expansion formula 
*22j plays an essential role to obtain Z n . We update 
11000 trajectories including 3000 thermalization ones. 
The measurement is performed at each 10 (20) trajctries 
for 8 3 x 4(10 3 x 4). The Monte Carlo update is performed 
with the fermion measure at ji = 0; thus we avoid any 
sign problem due to the complex fermion determinant. 
However, an overlap problem still exists, especially at 
low temperatures. When evaluating Z ni the correspond- 
ing operators fluctuate rapidly. 

The LYZ diagrams of the lattice QCD above and 
just below the phase-transition temperature are shown 
in Figs El and [71 respectively. For Z n , we use data of 
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FIG. 6: LYZ diagram at T/T c = 1.08 (/3 = 1.9) on the 
lattice. Lattice size is 10 3 x 4 and 8 3 x 4 and iV Max = 28 
and 14, respectively. 
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FIG. 7: LYZ diagram at T/T c ~ 1.0 fj3 = 1.85) on the 
QCD lattice. Lattice size is 8 3 x 4. 

(Signal/Noise) > 2. Since the simulation poorly performs 
at lower temperatures, fewer Z n points are available in 
those regions. The condition Z n = for n = l,2 (mod 3) 
(which can be easily proven) is imposed on the lattice 
data for Z„, which guarantees 27r/3 translational invari- 
ance for Im(/i)/T. We evaluate the LYZs for Z(£b) = 

Dm z 3m£,Bi an d ma P tne result onto £ = £ B ' . 

Note that, because the Z n obtained in the confine- 
ment region suffer from significant noise, the LYZ dia- 
gram should be considered qualitative. However, despite 
this, distinctive differences between Figs|6]and[7jare ap- 
parent. At T ^ 1-2T C , a line of zero accumulation ap- 
pears at arg(£) = ±7r/3, which indicates a Roberge- Weiss 
(RW) phase transition. 

Roberge and Weiss discuss the regions of pure imag- 
inary chemical potential, and find that at Im(/i)/T = 
7r/3 + 2fc7r/3 (fc = 0,1), a phase transition occurs for 
T > Trw |23j| . By using the same lattice setup as that 
used in the present study, the quantity Trw was esti- 
mated to be around 1.1T C [24]. This phase transition 
exerts a clear effect on the LYZ diagram in Fig|5] where 
the temperature T ~ 1.2T C . However, no such effect ap- 



4 



' Ful 


+ 




-1 






/ \-2 






/ % M \ 




/ * \ 




/ * \ 




/ i f \ 










! * 




\ * 












\ " « / 






i m r i i 






-1 -0.5 0.5 






Re[le] 





FIG. 8: LYZ diagram (STAR at y/s = 200 GeV). "Full" 
result is calculated by using all Z n . is by dropping 
a pair of Z n of n = ±Y max . For "-2", we drop further 
two pair Z n . The inset is an enlarged view of the zeros 
near positive real axis, which gives a predicted region of 
QCD phase transition shown in FigfJ] 



pears in Fig where the temperature is close to T c . 

We next consider the LYZ diagrams on the basis of 
RHIC Star data. We construct the grand partition func- 
tion Eq.® for RHIC energies, yfs = 19.6, 27, 39, 62.4 
and 200 GeV, from the data of P n > 2, and use the 
cBK method to calculate the LYZ diagrams. The result 
for y/s —200 GeV is shown in FigJ5] To make the finite 
]V m ax effect apparent, we also calculate the LYZs by ne- 
glecting Z n for n = ±N max ("-1") and by neglecting Z n 
for n = ±N max and ±(Y max - 1) ("-2"). 

Although some zeros exist on the negative real axis, 
these do not form a line that clearly characterizes the 
RW transition, which suggests that the data correspond 
to temperatures below T^w ■ Remember that the net- 
proton number is not an exactly conserved quantity, but 



it nevertheless may carry information at the time when it 
was born at higher temperature. Therefore, these points 
may be a trail of the line correspoiding to the RW tran- 
sition. 

In the LYZ diagrams for RHIC, no zeros appear in 
the positive real-£ region. Possible explanations for this 
result are that (i) there is no phase transition, but the 
cross-over occurs at these temperatures, (ii) the systems 
are finite, and/or (iii) the -/V max values are too small. 
The size of the fire ball produced is comparable to the 
QCD scale, and thus, explanation (ii) is at least partially 
correct. To further explore the QCD phase transition, 
larger A^nax must be attained. 

Finally, we estimate where the Lee- Yang zeros would 
intercept or reach near to the positive real axis as the 
volume increases. These zones are indicated by double 
-headed arrows in the inset of Fig(5] In FigfJ] the corre- 
sponding regions are indicated by horizontal lines. If the 
multiplicities P n with larger n were to be measured in 
future experiments, the QCD phase transition could be 
pinpointed with more precision. The further study of the 
relation between baryon and proton number distribution 
will improve the analysis here[25[. 
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